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We revise the analysis of the bottomonium hyperfine splitting within the lattice nonrelativistic 
QCD. The Wilson coefficients of the radiatively improved lattice action are evaluated by a semian- 
alytic approach based on the asymptotic expansion about the continuum limit. The nonrelativistic 
renormalization group is used to estimate the high-order radiative corrections. Our result for the 
IS" hyperfine splitting is M'f(is) — = 52.9 ± 5.5 MeV. It reconciles the predictions of the 

continuum and lattice QCD and is in very good agreement with the most accurate experimental 
measurement by the Belle collaboration. 

PACS numbers: 12.38.Gc, 12.38.Bx, 14.40.Pq, 14.65.Fy 


The bottomonium hyperfine splitting defined by the 
mass difference i?hfs = ATx(is) — has been a 

subject of much controversy since the first observation 
of the spin-singlet rji, state in radiative decays of the 
T(35') mesons by the Babar collaboration [l[. The mea¬ 
sured value of the hyperfine flitting overshot the predic¬ 
tions of perturbative QCD [2| by almost a factor of two, 
well beyond the experimental and theoretical uncertainty 
bands, see Table HI Further experimental studies [J-Q 
were consistent with the initial measurement, while the 
Belle collaboration reported a significantly lower value of 
the splitting with higher experimental precision Q. On 
the theory side the most accurate estimates of the hyper¬ 
fine splitting are obtained from lattice simulations within 
the effective theory of nonrelativistic QCD (NRQCD). 
This method is entirely based on first principles, al¬ 
lows for simultaneous treatment of dynamical heavy and 
light quarks and gives a systematic account of the long¬ 
distance nonperturbative effects of the strong interaction. 
The first analysis W with fully incorporated one-loop ra¬ 
diative corrections [8| favored the larger value of the split¬ 
ting fl. The most recent analysis includes the leading 
relativistic corrections an d g ives a lower value, which is 
close to the PDG average [lOj but nevertheless not consis¬ 
tent with Ref. 0 . This might indicate a serious failure of 
perturbative QCD in the description of the bottomonium 
ground state in clear conflict with the general concept of 
the heavy quarkonium dynamics. Thus the current ex¬ 
perimental and theoretical status of the bottomonium 
hyperfine splitting remains ambiguous and sets up one of 
the most interesting open problems in the QCD theory of 
hadrons, which yet inspired a discussion about possible 
signal of physics beyond the standard model [HI- 

In this article we revise the analysis of the radiative 
corrections to the lattice NRQCD action. We develop 
a semianalytical approach based on the asymptotic ex¬ 
pansion about the continuum limit [l^ . which provides 
a very powerful tool for the radiative improvement of 
lattice NRQCD. Our result for the one-loop Wilson co¬ 
efficient of the effective spin-dependent four-quark inter¬ 


action significantly differs from the result of the previ- 
ous calculation Q used in the subsequent analyses [3, , 

which leads to a sizable reduction of the lattice NRQCD 
prediction for the hyperfine splitting. We give an esti¬ 
mate of the higher order radiative corrections by evalu¬ 
ating the two-loop double-logarithmic terms within the 
nonrelativistic renormalization group approach 2, 141. 
The main result of this paper is a new theoretical value 
for bottomonium hyperfine splitting, Eq. (mi. 

The idea of the NRQCD approach [Hill is to sep¬ 
arate the hard modes, which require a fully relativistic 
analysis, from the nonrelativistic soft modes. The dy¬ 
namics of the soft modes is governed by the effective non¬ 
relativistic action given by a series in heavy quark veloc¬ 
ity V, while the contribution of the hard modes is encoded 
in the corresponding Wilson coefficients. The nonrela¬ 
tivistic action can be applied in a systematic perturba¬ 
tive analysis of the heavy quarkonium spectrum Yf Si- 
At the same time the action may be used for lattice sim¬ 
ulations of the heavy quarkonium states 0, im. The 
latter approach gives full control over nonperturbative 
long-distance effects and can be used for the description 
of excited states where perturbative QCD is not applica¬ 
ble. 

The hyperfine splitting i.e. the splitting between the 
spin-singlet and spin-triplet states is generated by the 
spin-dependent part of the NRQCD Lagrangian. To or¬ 
der Olv"^) it reads (see e.g. 22, 23|) 


Ca = + (V" ^ Xc) + (T'tpxlfTXc 

Z/TTtq TTLq 

( 1 ) 

where B is the chromomagnetic field, rriq and as are the 
heavy quark mass and the strong coupling constant, the 
SU{Nc) color group factor is Cp = — l)/{2Nc), 4> 

(Xc) are the nonrelativistic Pauli spinors of quark (anti¬ 
quark) held, and we have projected the four-quark inter¬ 
action on the color-singlet state. The Wilson coefficients 
cf and logarithmically depend on the factorization 
scale pif which separates the hard and the soft momen- 
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Experiment 

Babar, T(3S') decays[]J 

71.4l2;5(stat) ±2.7(syst) 

Babar, T(2S') decays [3] 

66 . 1/4 ® (stat) ± 2 . 0 (syst) 

Belle, hb{lP) decays 

57.9 ± 2.3 

PDG average [1^ 

62.3 ±3.2 

Theory 

NRQCD, NLL 

41 ± ll(th)/g((5as) 

Lattice NRQCD Oiv'^) [7] 

70 ±9 

Lattice NRQCD C>(a®) 

62.8 ±6.7 

Lattice QCD [3^ 

54.0 ± 12.4//® 

Lattice NRQCD, this work 

52.9 ±5.5 


TABLE I. Results of high-precision experimental and theoret¬ 
ical determinations of the bottomonium hyperfine splitting in 
MeV. 


turn contributions. This dependence can be predicted 
to all orders of perturbation theory by renormalization 
group methods. In lattice NRQCD the natural factoriza¬ 
tion scale is given by the inverse lattice spacing a. The 
radiative improvement of the action is therefore manda¬ 
tory for the correct continuum limit. 

The coefficient cf parametrizes the quark anomalous 
chromomagnetic moment. It can be determined nonper- 
turbatively by matching the lattice result for particular 
splittings to the physical bottomonium spectrum (?. 241. 
The perturbative evaluation of the one-loop correction 
to Cf [ 1 ] is in good agreement with the nonperturbative 
result. The Wilson coefficient of the effective four-quark 
interaction however can only be obtained perturbatively. 
It vanishes in the Born approximation and is determined 
by matching the one-particle irreducible quark-antiquark 
scattering amplitudes in QCD and NRQCD; see Fig. [T] 
The matching does not depend on the choice of soft kine- 
matical variables and 

becomes particulary simple when the amplitude is 
computed at the quark-antiquark threshold and vanish¬ 
ing momentum transfer. In this case the one-loop full 
QCD amplitude is 


m: 


QCD 

IPI 


Cpoii 


I - 


Ca 

2 

2'ITmq 

3A 


log -b (ln2 - l)rF 
Cf a-ipxlo-Xc, 


( 2 ) 


where Ca = Nc, Tf = 1/2, and we introduced a small 
auxiliary gluon mass A to regulate the infrared diver¬ 
gence. The power enhanced 1/A term corresponds to the 
Coulomb singularity of the threshold amplitude, while 
the term proportional to Tp is due to the two-gluon an¬ 
nihilation of the quark-antiquark pair. 

On the other hand the lattice NRQCD result for the 



FIG. 1. One-loop Feynman diagrams contributing to the one- 
particle irreducible quark-antiquark scattering amplitude in 
QCD [(a)-(d)] and NRQCD [(c)-(f)]. 


one-loop amplitude can be written as follows 


a^NRQCD _ CfoI 

iWipi — 

rriq 
27r rriq 


Cf 


^5+iln(aA)^ Ca 
ct'4}xIctxc + O ( a ^) , 


( 3 ) 


where the nonlogarithmic non-Abelian term 5 depends on 
a particular realization of the lattice action. To match 
Eqs. ([2]) and ([3]) we add to the NRQCD Lagrangian the 
four-quark operator with coefficient 


djfj — 


i5 -l- — T j Ca + (la 2 — 1 )Tf -I-Cf 


( 4 ) 


where L = In(amq). The main problem is therefore in 
determination of the coefficient 5. The asymptotic ex¬ 
pansion of the lattice loop integrals about the continuum 
limit 12| can in principle be used to get this coefficient 
in a closed analytic form. Since the heavy quark mass is 
not a dynamical scale in NRQCD, the parameter of the 
expansion in our case is a\. The idea of the method is 
to split the integration over the virtual momentum I into 
the contributions of the hard region with Z ^ 1/a and the 
soft region with Z ~ A. In the hard region the integrand is 
expanded in a A and A/Z and reduces to the lattice tadpole 
integrals. In the soft region the integrand is expanded in 
al and aX and reduces to the continuum NRQCD Feyn¬ 
man integrals. As a result of the scale separation the 
hard (soft) contribution in general has spurious infrared 
(ultraviolet) logarithmic divergences and has to be regu¬ 
lated. In the total result for a given lattice loop integral 
the dependence on the regulator cancels out leaving the 
asymptotic series in aX which includes the logarithmic 
terms [cf. Eq. ([3])]. We emphasize that the expansion 
about the continuum limit is a formal tool to get the 
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NRQCD loop integrals as series in a and facilitate the 
matching procedure, while the lattice NRQCD is a valid 
nonrelativistic effective theory only for a ^ ^Irriq. Note 
that Eq. ([31) has only a logarithmic singularity in a. In 
higher orders of the NRQCD expansion in l/iriq the am¬ 
plitude includes also the terms with a negative power of 
a. Such l/(amg)" terms are more singular in the formal 
continuum limit but are power suppressed with respect 
to Eq. dH) in the region where lattice NRQCD is applied. 

Let us consider first a “naive” lattice action with no 
improvement for gluonic and heavy quark fields (see, e.g. 


26|1. The gluonic field tensor of the NRQCD chro- 
momagnetic interaction in the naive action is expressed 
through the commutator of the left-right symmetrized 
covariant lattice derivatives. In this case we obtain 

^naive = -1 + 28 n'^b 2 - = 0.288972 ..., (5) 


where the irrational constants 62 = 0.02401318..., 63 = 
0.00158857... parametrize the lattice tadpole integrals 
and can be computed with arbitrary precision [l^ . We 
however need the above coefficient for the improved lat¬ 
tice action which is used in real simulations. Analytic 
calculation with an improved action is not optimal since 
the Feynman rules in this case become extremely cum¬ 
bersome. We bypass this problem by using a semiana- 
lytic approach. Indeed the difference between the Wil¬ 
son coefficients for the improved and naive lattice ac¬ 
tions Ai5 remains finite in the limit A —>■ 0 and can be 
directly obtained by numerical evaluation of the corre¬ 
sponding lattice loop integrals with sufficiently small A 
(a finite infrared regulator is necessary for the stabil¬ 
ity of numerical integration). For the numerical imple¬ 
mentation of the improved lattice action Feynman rules 
we use HiPPy/HPsrc code However in contrast 

to the standard implementation the color space reduc¬ 
tion is performed analytically with the help of the pro¬ 
gram COLOR [ 2 ^ before the numerical integration is 
done by the CUBA integrator library [ 2 ^. This greatly 
reduces the runtime and allows for a separate treat¬ 
ment of the contributions of independent color group 
structures which have different infrared properties, c/. 
Eq. ([3|). The whole process of the calculation is fully au¬ 
tomated. In the case of the HPQCD action Q we get 
A(5 = —0.1444(28) corresponding to 


<5 = 0.1446(28). (6) 


Note that since do- = 0 in the Born approximation, we 
have to perform neither the strong coupling constant 
renormalization nor the lattice tadpole improvement. We 
made a few cross-checks of the calculation. For the naive 
action the numerical evaluation agrees with the gauge- 
invariant analytic result of the asymptotic expansion for 
small values of A. The logarithmic part of is in agree¬ 
ment with the renormalization group analysis. The non¬ 
relativistic renormalization group predicts the all-order 





FIG. 2. The result of the lattice simulation of the bottomo- 
nium hyperfine splitting with 0{v^) action [J and 0{v^) ac¬ 
tion [^. The error bars are explained in the text. The solid 
lines correspond to the central values of the constrained fit. 


dependence of the Wilson coefficients on (Hill. In 
the leading logarithmic approximation they read 


= 


Ca 


/3o - 2Ca 




„LL _ -Ca 
C p — z , 


(7) 


where /3o = llC^/d — 4Tpni/3 is the one-loop QCD 
fd function, n/ = 4 is the number of light flavors, and 
z = {as{gif)/as{mq)Y/^°. In lattice NRQCD the fac¬ 
torization scale should be identified with inverse lattice 
spacing /r/ ^ 1/a. By reexpanding the leading logarith¬ 
mic result we obtain 


jLL 


Ois Ca 


L- 




(^0 


Ca)Ca ^2 


Ti = 1 - 


Cp 


tts ^ (/3o + Ca) Ca j 2 

- 1j 


m 


( 8 ) 


in agreement with Eq. (SD- 

Let us now compare our result with the previous calcu¬ 
lation Q . In this paper a different basis of the four-quark 
operators is used and the Wilson coefficient da jots should 
be identified with the linear combination |(di — ^ 2 ) (see 
( 3 ^ for the consistent analytical expressions). We find 
that the nonannihilation constant term of the QCD am¬ 
plitude in Q is smaller than the one in Eq. m by a factor 
of 3. The comparison of the NRQCD part of the result is 
not straightforward since in [8|, l30| it has been evaluated 
numerically for three different lattice spacings keeping 
the full dependence on rriq. This dependence includes 
power suppressed terms as well as the linear term from 
the lattice cutoff of the Coulomb pinch contribution. For 
the lattice spacing used in real simulations a ^ l/{vmq) 
the power suppressed terms are of the same magnitude 
as the generic one-loop relativistic corrections and are 
beyond the accuracy of our analysis. On the contrary 
the lattice artifacts are numerically significant. The lin¬ 
ear term associated with the Coulomb pinch can be esti¬ 
mated by cutting the corresponding continuum NRQCD 
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0(v'^) action 

0{v^) action 

Discretization error 

2.6 

3.1 

Relativistic corrections 

6.0 

1.8 

Radiative corrections 

4.8 

4.3 

-E^hfs 

57.5 

51.5 


TABLE II. The central value and the error budget for the 
lattice NRQCD determination of the bottomonium hyperfine 
splitting with 0(v'^) action [J and 0(v^) lattice action in 
MeV. 


one-loop integral at the scale tt/o. This gives an addi¬ 
tional contribution to Eq. o, 


8 C 

— u - amn 

3 TT ^ 


—0.94 ttsamq , 


(9) 


where the factor v = 0.831... adjusts the analytical re¬ 
sult for the integral over spherical momentum domain 
to the integral over the Brillouin zone. The numerical 
result of Ref. Q suggests a significantly larger negative 
coefficient of about —1.8. Moreover in the threshold re¬ 
gion the multiple Coulomb gluon exchange contributions 
are not parametrically suppressed and the modification 
of the Coulomb bound state dynamics on the finite lat¬ 
tice is not accounted for by the one-loop analysis. It may 
change the numerical coefficient of the linear term and 
generates all-order contributions in aganiq. This means 
that (i) the one-loop matching does not remove the lin¬ 
ear lattice artifact at 0(as) and (ii) one cannot use the 
finite lattice spacing a ^ l/iasTUq) as a Wilsonian cut¬ 
off for NRQCD as it was done in [§, |3^. Thus all the 
lattice artifacts should be removed nonperturbatively W 
numerical extrapolation of the lattice data to a = 0 [7|,[9|. 

Now we are in a position to apply our result to the 
analysis of the hyperfine splitting. The contribution of 
the four-quark interaction to Ejus reads 


— d(j 


^CpOis 

ml 


\m\^ 


( 10 ) 


where is the wave function of the quarkonium 

ground state at the origin. Equation should be 
added to the result of the lattice simulation with the 
one-loop Wilson coefficient cp and no four-quark contri¬ 
bution included. Such a result is available for the 
action Q and for the 0{v^) action [^. For the numerical 
analysis of Eq. (|10|) we use the nonperturbative lattice 
result for [7|. To make our analysis self consistent 
we adopt the value of the bottom quark mass to;, and 
the value of the strong coupling constant ay renormal¬ 
ized in the static potential scheme at the scale n/a from 
Ref. 0- The numerical result for the hyperfine splitting 
is presented in Fig. [2] as a function of for three dif¬ 
ferent lattice spacings and two different lattice actions. 
The error bars of each point include the statistical er¬ 


ror and the uncertainty in the value of the lattice spac¬ 
ing from 0,0 as well as the high-order a-dependent ra¬ 
diative corrections that are estimated by the size of the 
double-logarithmic two-loop terms in Eq. ([8]). The use of 
relatively large values of the lattice spacing a l/{vmb) 
ensures the suppression of the unphysical l/(aTO{,)” con¬ 
tributions, which become important at a ^ l/mb 00. 
At the same time it results in sizable lattice artifacts, 
which cannot be removed by finite order matching due 
to the all-order character of the Coulomb binding effects. 
To minimize this effect the result is numerically extrapo¬ 
lated to o = 0 00. The extrapolation below a ^ 1 /mb 
in this case is justified since the numerical effect of the 
l/(aTOh)"' terms on the data points is small. To perform 
the extrapolation we use a constrained fit of the data 
points [3l| by a polynomial in a with vanishing linear 
term. The inclusion of the linear and l/(amf,)” terms in 
the fit is discussed below. To estimate the coefficients of 
the higher order terms in the lattice spacing we represent 
the result of the fit as 1 -I- (Aa)^ -I- 0{a^), where A is the 
mass scale characterizing the approach of the lattice ap¬ 
proximation to the continuum limit. The priors for the 
coefficients of the a" terms with n > 2 in the constrained 
fit are then given by the intervals ±A". Numerically we 
get A « 360 MeV for the 0{v‘^) and A m 790 MeV for the 
0{v^) case. Because of a slower approach to the contin¬ 
uum limit the extrapolation error for O(v^) action turns 
out to be larger. This may be related to the fact that the 
contribution of the operators of higher dimension is more 
sensitive to the ultraviolet momentum region. Therefore 
the currently unknown O(asV^) matching corrections in 
this approximation can be substantial. We checked that 
the inclusion of the 1/a" terms with the priors 
into the constrained fit changes the result within the ex¬ 
trapolation error intervals. 

In general the Coulomb binding effects give rise to a 
linear dependence of the lattice data on a which can be 
roughly estimated by the one-loop result . A more 
refined estimate can be obtained by including the linear 
term ciagamq into the fit of the lattice data. For the 
prior |c/| < 1 the constrained fit gives c/ ~ —0.25 for 
both actions, which is two times smaller than the one- 
loop estimate ci w —0.5 corresponding to Eq. Q. At the 
same time the extracted value of the hyperfine splitting 
is increased within the extrapolation error interval by 
approximately 2.5 MeV. 

The total error budget of our estimate is given in Ta¬ 
ble El Besides the discretization errors discussed above 
it includes the uncertainty due to high-order relativistic 
and radiative corrections. For a conservative estimate of 
the radiative corrections we take the value of the double- 
logarithmic two-loop terms at the soft factorization scale 
/j,/ aemb dictated by the bound state dynamics. In 
Table El this uncertainty is combined with the numerical 
error in the one-loop coefficient cp 0 . Our estimate of 
the relativistic corrections for the 0(v‘^) action is based 
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on the difference between the 0{v^) and 0{v^) results in 
the continuum limit. For the 0{v^) action we multiply 
this uncertainty by ag evaluated at the soft renormaliza¬ 
tion scale to take into account the previously discussed 
missing matching corrections. The larger discretization 
uncertainty balances the smaller relativistic corrections 
in the 0{v^) case and both actions provide comparable 
total errors. Since the structure of the relativistic cor¬ 
rections and the behavior of the results at finite lattice 
spacing are significantly different for the two actions, we 
consider the corresponding uncertainties as uncorrelated 
and take the weighted average of the results as the best 
estimate. At the same time the uncertainty due to the 
high-order purely radiative corrections is treated as cor¬ 
related between the two actions. Our final result for the 
hyperfine splitting reads 

Ahfs = 52.9 ± 5.5 MeV. (11) 

We now can compare our estimates to the available the¬ 
oretical and experimental results in Table HI Our result 
for both 0{v*) and 0{v^) actions fTable 111)) are below 
the ones of the previous lattice NRQCD analysis [ 3 , 0 
by approximately 12 MeV. About 5 MeV of the differ¬ 
ence is due to the error in the one-loop QCD amplitude 
calculation [^. The remaining discrepancy is related to 
the different procedure of extrapolation to a = 0. The 
analysis 0,0 implies that the one-loop matching 0 re¬ 
moves the linear artifact form the lattice data. However, 
as it was pointed out above, the one-loop calculation can 
only be used for a rough estimate of the linear term due 
to the all-order character of the Coulomb binding effects. 
We therefore determine the corresponding coefficient by a 
constrained fit of the lattice data with the prior set by the 
one-loop result. Moreover the numerical result 0 sug¬ 
gests a significantly larger value of the linear term than 
what follows from our analytic calculation and from the 
fit of the lattice data, which leads to a sizable difference 
of the extrapolation results. 

With the new value of the four-quark Wilson coeffi¬ 
cient the lattice NRQCD prediction m agrees within 
the error bars with the next-to-leading logarithmic (NLL) 
perturbative QCD result 0 . Its central value practically 
coincides with that of the full lattice QCD simulation 
[^ . though the uncertainty of the latter is significantly 
larger. This may indicate that the matching of the lat¬ 
tice NRQCD to full QCD is now done properly. On the 
experimental side our result strongly favors the value ob¬ 
tained by the Belle collaboration, which has the lowest 
reported uncertainty. Thus we have reconciled the theo¬ 
retical predictions of the lattice and continuum QCD as 
well as the most accurate experimental data. 
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